In a previous notebook, I have shown that the . Here, I would like to use the existing database of stars' positions and display them as a triangulation

Let's first initialize the notebook:

In [12]:

import numpy as np
np.set_printoptions(precision=6, suppress=True)
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
phi = (np.sqrt(5)+1)/2
fig_width = 10
figsize = (fig_width, fig_width/phi)

importing data¶

Importing all stars' position is as simple as invoking the HYG database:

In [118]:

import pandas as pd
url = "https://github.com/astronexus/HYG-Database/raw/master/hygdata_v3.csv"
space = pd.read_csv(url)

In [119]:

print(f'Columns in the database: {space.columns=}')

Columns in the database: space.columns=Index(['id', 'hip', 'hd', 'hr', 'gl', 'bf', 'proper', 'ra', 'dec', 'dist',
       'pmra', 'pmdec', 'rv', 'mag', 'absmag', 'spect', 'ci', 'x', 'y', 'z',
       'vx', 'vy', 'vz', 'rarad', 'decrad', 'pmrarad', 'pmdecrad', 'bayer',
       'flam', 'con', 'comp', 'comp_primary', 'base', 'lum', 'var', 'var_min',
       'var_max'],
      dtype='object')

In [120]:

print(f'Number of stars in the catalog = {len(space)=}')

Number of stars in the catalog = len(space)=119614

extracting ra, dec and mag¶

For which we may extract what interests us: position (right ascension and declination) and visual magnitude:

In [199]:

space_pos = space[['ra', 'dec', 'mag']]
space_pos.is_copy = False
space_pos

Out[199]:

	ra	dec	mag
0	0.000000	0.000000	-26.70
1	0.000060	1.089009	9.10
2	0.000283	-19.498840	9.27
3	0.000335	38.859279	6.61
4	0.000569	-51.893546	8.06
...	...	...	...
119609	23.963895	38.629391	12.64
119610	23.996567	47.762093	16.10
119611	23.996218	-44.067905	12.82
119612	23.997386	-34.111986	12.80
119613	0.036059	-43.165974	13.05

119614 rows × 3 columns

First, right ascension is "the angular distance of a particular point measured eastward along the celestial equator from the Sun at the March equinox to the (hour circle of the) point in question above the earth." It is given in hours:

In [200]:

ra_min, ra_max = 0, 24
print(f"RA: {space_pos['ra'].min()=}, {space_pos['ra'].max()=}")

RA: space_pos['ra'].min()=0.0, space_pos['ra'].max()=23.998594

Let's convert this in visual angle, that is in an azimuth:

In [201]:

az_min, az_max = 0, 360
def ra2az(ra):
    return az_max - ra / ra_max * az_max
space_pos['az'] = ra2az(space_pos['ra'])
print(f"AZ: {space_pos['az'].min()=}, {space_pos['az'].max()=}")

AZ: space_pos['az'].min()=0.021090000000015152, space_pos['az'].max()=360.0

/var/folders/3p/m0g52j9j69z3gj8ktpgg1dm00000gn/T/ipykernel_6817/1446867934.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  space_pos['az'] = ra2az(space_pos['ra'])

Then, declination is "comparable to geographic latitude, projected onto the celestial sphere" and is given here in degrees:

In [202]:

dec_min, dec_max = -90, 90
print(f"DEC: {space_pos['dec'].min()=}, {space_pos['dec'].max()=}")

DEC: space_pos['dec'].min()=-89.782428, space_pos['dec'].max()=89.569427

The magnitude varies a lot (the less the value, the more it is visible - "The apparent magnitudes of known objects range from the Sun at −26.7 to objects in deep Hubble Space Telescope images of magnitude +31.5"):

In [203]:

print(f"mag: {space_pos['mag'].min()=}, {space_pos['mag'].max()=}")

mag: space_pos['mag'].min()=-26.7, space_pos['mag'].max()=21.0

Let's normalize the lower bound:

In [204]:

space_pos['mag'] = space_pos['mag'] - space_pos['mag'].min()
print(f"mag: {space_pos['mag'].min()=}, {space_pos['mag'].max()=}")

mag: space_pos['mag'].min()=0.0, space_pos['mag'].max()=47.7

/var/folders/3p/m0g52j9j69z3gj8ktpgg1dm00000gn/T/ipykernel_6817/2801712572.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  space_pos['mag'] = space_pos['mag'] - space_pos['mag'].min()

Let's define a threshold using the quantile function:

In [205]:

print(f"{space_pos['mag'].quantile(q=0.01)=}")

space_pos['mag'].quantile(q=0.01)=31.45

Which leaves only a limited number of stars

In [206]:

space_bright = space_pos[space_pos['mag']<space_pos['mag'].quantile(q=0.05)]

In [207]:

#space_bright = space_pos[space_pos['mag']<6]

In [208]:

print(f'Number of bright stars = {len(space_bright)=}')

Number of bright stars = len(space_bright)=5955

In [209]:

print(f"MAG: {space_bright['mag'].min()=}, {space_bright['mag'].max()=}")

MAG: space_bright['mag'].min()=0.0, space_bright['mag'].max()=32.85

scatter plots¶

From these elements, we may plot the stars on these coordinates:

In [210]:

fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.scatter(space_bright['az'], space_bright['dec'], s=.1 * (space_bright['mag'].max()-space_bright['mag']))
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(az_min, az_max)
ax.set_ylim(dec_min, dec_max);

do you see the Milky way?

Orion¶

We may want to focus on the Orion constellation):

In [237]:

orion_az_max, orion_az_min = ra2az(4.8), ra2az(6.25)
orion_dec_min, orion_dec_max = -13.0, 22.
if True:
    space_orion = space_bright[(orion_az_min < space_bright['az']) * (space_bright['az'] < orion_az_max) *
                                (orion_dec_min < space_bright['dec']) * (space_bright['dec'] < orion_dec_max)]
else:
    space_orion = space_pos[(orion_az_min < space_pos['az']) * (space_pos['az'] < orion_az_max) *
                            (orion_dec_min < space_pos['dec']) * (space_pos['dec'] < orion_dec_max)]

In [238]:

print(f'Number of bright stars in orion = {len(space_orion)=}')

Number of bright stars in orion = len(space_orion)=192

In [239]:

print(f"MAG: {space_orion['mag'].min()=}, {space_orion['mag'].max()=}")

MAG: space_orion['mag'].min()=26.88, space_orion['mag'].max()=32.85

In [240]:

fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.scatter(space_orion['az'], space_orion['dec'], s= 10 * (space_orion['mag'].max()-space_orion['mag']))
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(orion_az_min, orion_az_max)
ax.set_ylim(orion_dec_min, orion_dec_max);

Which is to be compared to this image:

In [242]:

from  IPython.display import display, SVG, Image
display(SVG('https://upload.wikimedia.org/wikipedia/commons/f/ff/Orion_IAU.svg'))

Triangulation¶

Following work with Etienne Rey, we have developped visualization based on triangulation of an optimal packing of point clouds.

In [243]:

display(Image('https://laurentperrinet.github.io/post/2021-10-04_interstices/featured.jpg'))

The idea here is to apply such a visualization to the stars position. IT is based on

In [245]:

import matplotlib.tri as tri
triang = tri.Triangulation(space_bright['az'], space_bright['dec'])

In [247]:

fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.triplot(triang, lw=0.9, color='k')
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(az_min*.9, az_max*.9)
ax.set_ylim(dec_min*.9, dec_max*.9);

In [249]:

fig, ax = plt.subplots(figsize=(fig_width, fig_width))
ax.triplot(triang, lw=0.9, color='k')
ax.scatter(space_orion['az'], space_orion['dec'], s= 40 * (space_orion['mag'].max()-space_orion['mag']))
ax.set_aspect('equal')
ax.set_xlabel('Azimuth')
ax.set_ylabel('Declination')
ax.set_xlim(orion_az_min, orion_az_max)
ax.set_ylim(orion_dec_min, orion_dec_max);

some book keeping for the notebook¶

In [ ]:

%load_ext watermark
%watermark -i -h -m -v -p numpy,matplotlib,scipy,pillow,imageio  -r -g -b